lightfield = zeros(16,16,400,700,3);
lightfield(i,j,:,:,c) = lf_img(i:16:6400,j:16:11200,c);
t = tiledlayout(16, 16, 'TileSpacing', 'none', 'Padding', 'compact');
img = uint8(squeeze(lightfield(i,j,:,:,:)));
img = squeeze(lightfield(i,j,:,:,:));
t = tiledlayout(1,1, 'TileSpacing', 'none', 'Padding', 'compact');
imshow(uint8(imgs/(16*16)))
d_list = {0.25, 0, -0.25, -0.5, -0.75, -1, -1.2, -1.5};
img = squeeze(lightfield(i,j,:,:,:));
s = round(S-d_list{d}*u(i)); s(s>400) = 400; s(s<1) = 1;
t = round(T+d_list{d}*v(j)); t(t>700) = 700; t(t<1) = 1;
imgs = imgs + img/(16*16);
tmp = rgb2xyz(fstack{i});
fstack_lum{i} = tmp(:,:,2);
%2 Create low frequency component
tmp = im2double(imgaussfilt(fstack_lum{i},sigma1));
%3 Compute high frequency component
fstack_high{i} = fstack_lum{i} - fstack_low{i};
%4 Compute the (sharpness) weight
tmp = im2double(imgaussfilt((fstack_high{i}).^2,sigma2));
t = tiledlayout(1,1, 'TileSpacing', 'none', 'Padding', 'compact');
nom = nom + wsharp{i}.*im2double(fstack{i});
denom = denom + wsharp{i};
focused_img = nom./denom;
nom = nom + wsharp{i}.*(d_list{i}+1.5);
denom = denom + wsharp{i};
imshow((depthmap-min(depthmap,[],'all'))/(max(depthmap,[],'all')-min(depthmap,[],'all')));